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Abstract 

We present three-dimensional models of accretion disks in U Gem - 
like systems and calculate their Doppler tomograms. The tomograms are 
based on two different assumptions concerning the origin of line emission 
from the disk. The assumption of lines originating due to irradiation of the 
surface layer of the disk by the central source leads to a better agreement 
with observations. We argue that fully three-dimensional modelling is 
necessary to properly interpret the observed tomograms. 

1 Introduction 

With their typical dimensions of less than 10~* arcseconds, Dwarf Nova 
(DN) disks are far too small to be resolved directly. Ifowever, an indi- 
rect observational insight into their structure became possible already in 
mid-eighties, when a powerful technique of Doppler tomography was in- 
troduced (for a recent review see Marsh 2001). Nowadays tomographic 
observations are a standard, but the interpretation of Doppler tomograms 
in terms of theoretical models of DN disks is often problematic. This is 
particularly well visible in the case of spiral waves, which have been sug- 
gested as one of the agents responsible for the angular momentum transfer 
through the disk (for a recent review see BofBn 2001). The excellent dis- 
cussion of the subject can be found in a recent paper by Smak (2001), 
and there is no need to repeat it here. 

The main problem associated with the interpretation of Doppler tomo- 
grams is related to the nature of the data derived from the theory. Model 
calculations yield spiral patterns in the distribution of disk surface density, 
while patterns in observational tomograms are related to the distribution 
of the emissivity in specific spectral lines. While comparing these two 
sets of data one usually makes an implicit assumption that emissivity is 
proportional to surface density, which certainly is an oversimplification. 

A more sophisticated approach was presented by Steeghs and Stehle 
(1999), who based their Doppler tomograms on emission line profiles cal- 
culated from 2D disk models. For the origin of the lines they adopted a 
purely thermal model with local Planckian source function. Such model 
was criticized by Smak (2001), who pointed out that it requires factor 
of 100 overabundance of helium in order to reproduce the observed in- 
tensities, while the relative brightness of the features it produces in the 
tomograms is incompatible with observations. 

Smak himself proposed to base the theoretical tomograms on the dis- 
tribution of velocity divergence. He argued that at a given radial distance 
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from the disk center the compression regions with Vv < would be dis- 
tinguished by a higher-than-average disk thickness (because their density 
and temperature would be higher). As a result, the surface layer of the 
disk would be better exposed to the irradiating flux from the white dwarf 
and the boundary layer, and a local enhancement in line emission would 
be observed. Based on the three-body model of gas flow in a close binary 
he identified regions of maximum compression in the orbital plane and 
showed that they well reproduced shape, location and relative intensities 
of the arch-like structures observed in Doppler tomograms of DN disks. 

Motivated by his paper, we obtained two- and three-dimensional hy- 
drodynamical models of DN disks and calculated their tomograms. The 
details of the modelling procedure are given in Section |21 The results are 
presented in Section |3| and discussed in Section 0| 

2 Numerical methods, input physics and 
initial conditions 

All models presented here were obtained with the help of the ZEUS-3D 
code (Clarke & Norman 1994, Clarke 1996). The original code was mod- 
ified to include conservative angular momentum transport (Kley 1998). 
Spherical coordinates (r, 6, 0) centered on the primary and corotating with 
the system were used. The grid was extending from to 27r in (j), from 
to 0.2n in 9, and from rin — 0.1a to rout ~ 0.5a in r, with a standing for 
the orbital separation. Grid spacing was uniform in {9, 0) and logarith- 
mic in r, resulting in zones of identical shape. After a few experiments we 
decided to limit the resolution to 100 x 20 x 64 zones in r, 9, <j) directions, 
respectively (test runs, with resolution increased by a factor of 2 in r, 9, 
consecutively, did not introduce any significant changes into the models) . 

A standard periodic boundary condition was imposed ai (j) — 2tt, and 
symmetry with respect to the orbital plane was assumed, implying a re- 
flecting boundary condition at ^ = 0. A free outflow was allowed for at 
Tin and rout. In the four innermost radial zones the radial velocity was 
reduced by 5% at every time-step, preventing the reflection of waves from 
the inner boundary of the grid. To check whether this damping procedure 
did not influence the structure of the disk or the shape of the spiral pat- 
tern, we calculated model B2 with rin moved to 0.045a (see Table Ql. It 
was found that shifting the inner grid boundary toward the white dwarf 
did not introduce any signiflcant changes in the model. 

The simulations did not include explicit viscosity (the von Neumann & 
Richtmyer and scalar linear artiflcial viscosities originally implemented in 
ZEUS were only used with coefficients Ci = 0.5 and C2 ~ 2.0, where Ci 
is responsible for the magnitude of the artificial viscous pressure, and C2 
is a shock-spreading parameter, see the definitions in Stone and Norman 
(1992)). The energy equation was not solved; a polytropic equation of 
state (p — Kp^) with 7 = 5/3 was employed instead. In a system of units 
in which gravitational constant, orbital separation, and primary's mass 
are all equal to 1, the value of k was set to 6500. To mimic U Gem - like 
systems, all models had the same mass ratio, q = 0.5. The stream flowing 
from the secondary through the Li point was not included. 

Every simulation consisted of three phases (relaxation, switch-on, and 
proper). The mass ratio was set to in the relaxation phase and to 0.5 in 
the proper phase, while in the switch-on phase it was linearly increasing in 
time. At the beginning of each simulation {t — 0) the grid was initialized 
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Table 1: List of models 



Model 


type 


^radial 


r-in 


remarks 


A 


2D 


100 


0.1a 




Bl 


3D 


100 


0.1a 




B2 


3D 


150 


0.045a 


outer 100 zones placed as in Bl 



with an exponential density distribution 

p{r,0) = max (^poe~'"' "jPminj, (1) 

where 

- ^""^ \ (2) 



<o' 



and 

.2 dp 



dp 



= l^Po • (3) 

e=o 

In our system of units the value of the midplane density, po, was set to 
10^*, eorresponding to ~ 2 ■ 10^*^ g cm^'^ in U Gem (with U Gem parame- 
ters taken from Groot 2001). The azimuthal velocity of the disk was given 
a purely Keplerian pattern, and the remaining two velocity components 
were set to 0. Since the content of the grid was not in hydrostatic equilib- 
rium, we allowed it to relax for ~ 3.9 orbital periods (Porb). Throughout 
the relaxation phase the model was strictly axisymmetric, so that it was 
possible to speed the computations up by reducing the number of angular 
grid points to 2. 

Because of extremely steep vertical gradients of density at the surface 
of the disk it was necessary to introduce a density limit. Every time step 
the grid was scanned for cells with p < pmin, and whenever such a cell 
was found p was reset to pmin. We wanted pmin to be small enough to 
minimize side-effects caused by the newly-added matter falling onto the 
disk, and, simultaneously, largo enough to avoid excessive computational 
slowdowns duo to formation of strong shocks in the rarefied medium above 
the disk. After a few experiments pmin was set to 10~^^ for all models. At 
the end of the relaxation phase the midplane density of the disk increased, 
reaching up to 4 ■ 10~* (corresponding to ~ 8 • 10~* gcm~^ in U Gem). 

At the beginning of the switch-on phase the relaxed model was mapped 
onto the standard grid, and the secondary's gravity was "switched on". 
The final value of q was achieved at t = 5.8 Porb. At the end of the switch 
on phase the total mass contained in the grid, and scaled to U Gem, 
A/disk, was equal ~ 10'^'* g). The proper phase with q = 0.5 lasted for 
another ~ 3.9 Porb so that the simulation was ended at t ~ 9.7 Porb. By 
that time shape and location of the disk edge stabilized, and a stationary 
spiral pattern developed in the disk. 



3 Results 

The tidal forces affect the relaxed disk in two ways. First, some material 
is stripped from its outer edge and driven out of the grid through the 
outer grid boundary. Second, angular momentum is removed from the 
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Figure 1: Left panel: distribution of disk density, p, in model A. Right panel: 
distribution of disk surface density, E = J p{z)dz, in model Bl. The white loop 
marks the Roche lobe. 



remaining material and transferred into the orbital momentum of the 
binary, causing the disk to shrink. In the three-body approximation, 
the radius of the disk cannot be larger than the radius of the largest 
non-intersecting orbit, r,™™. U Gem - like systems with q — 0.5 have 
rJJIg" ~ 0.3, and in fact at the end of the simulation the disk barely extends 
beyond r ~ 0.3. The rest of the gas originally located at 0.3 < r < 0.5 now 
resides in a ring-like density enhancement between r ~ 0.15 and r ~ 0.3. 

The ring is markedly elliptical, but, when averaged over the azimuthal 
angle, it shows a well-defined density maximum at r ~ 0.19, i.e. slightly 
beyond the circularization radius (rdrc = 0.16 for q = 0.5). The maximum 
density is factor of ~ 2.5 higher than the midplane density of the inner 
disk {r < 0.15). The ratio h/r, where h is the half-thickness of the disk, 
varies from ~ 0.1 in the inner disk to ~ 0.2 in the ring. The overall 
structure of the final model is reminiscent of the one expected at the early 
phase of the outburst, when the outer radius of the disk just begins to 
increase. However, because of the polytropic equation of state we employ, 
the interior of the disk is unrealistically hot (for a disk composed of pure 
hydrogen T grows from ~ lO'' K at the surface of the disk to ~ 6.7 • 10^ K 
at the density maximum). We find that both 2-D and 3-D calculations 
produce disks of nearly the same shape and extent (Fig. ^ . Below we 
shall argue, however, that fully three-dimensional models are needed to 
properly interpret the Doppler tomograms. 

Both 2-D and 3-D models develop spiral shocks shown in Figs.QandlJl 
At a first glance, location and inclination of the shocks do not significantly 
depend on the number of dimensions of the model. Significant difi^erences 
become visible when compression regions (Vu < 0) are compared: while 
model A has a clear two-armed pattern, three arms are present in mod- 
els Bl and B2. We speculate that the third arm may originate due to 
tidal forcing of the disk matter in the direction perpendicular to the or- 
bital plane (an effect entirely absent in 2D). Thus, to the long dispute 
on whether spiral shocks can exist in three dimensions, or rather their 
existence is limited to the two-dimensional world, we add a vote in favor 
of the first possibility. In 3D the shocks are definitely there, but their 
pattern is different than in 2D. Obviously, the validity of this conclusion 
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Figure 2: Left panel: distribution of the comprcssional power per unit volume, 
—min{p\Iv, 0), in model A. Right panel: distribution of the compressional power 
per unit surface, — J mm(pV'u, 0)dz, in model Bl. 



is limited to hot polytropic disks with 0.1 < h/r < 0.2. 

Following the approach of Smak (2001), we obtained Doppler tomo- 
grams of the compressional power due to tidal forces, pdV .The distribu- 
tion of brightness on the {vx,Vy) plane was calculated from the integral 

Lpdv{vx,Vy) = ~ j rmn{pVu,Q)B{ux,Va:,S^^)B{uy,Vy,Syy)dV, (4) 
Jv 

where {ux,Uy) are velocity components of the volume element AV , and 
the so-called boxcar function B is defined as 

^(-'"''^) = { i l"erwi^e^ 

The resolution parameter 5, related to the resolution of images (400 x 400 
pixels), was set to 6t)orb/400. 

In the three-body approximation of Smak (2001), the inner disk is 
essentially Keplerian, and its Doppler tomogram cannot show any non- 
axisymmetric structures in the high-velocity range. However, in a more 
realistic polytropic disk the spiral shocks excited at its outer edge prop- 
agate deeply into the high-velocity regions. As a result, our tomograms 
(Fig. 01 show extended spirals instead of two crescent-shaped maxima re- 
ported by Smak (2001). Such spirals are not observed in real disks (Groot 
2001). To account for the limited resolution of the observational data we 
blurred the tomograms to such a degree that the width of the brightest 
segments of the spiral became comparable to the width of the observed 
features (~ 400 ~ 500kms~^; see Groot 2001). However, the spiral pat- 
tern extending up to velocities of ~ 1000 kms~^ was still clearly visible. 
Moreover, both location and relative intensity of the brightest segments 
of the spiral did not agree with those observed in U Gem. 

Should this disagreement be regarded as an argument against the pres- 
ence of spiral waves in CV disks? Certainly not. As we already indicated 
in the Introduction, the observed tomograms refer to the distribution of 
the emissivity in specific spectral lines rather than to the distribution 
of physical parameters directly obtainable from hydrodynamical simula- 
tions. The presently available hydrocodes are not sophisticated enough 
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Figure 3: Doppler tomograms of compressional power distributions shown in 
Fig. 13 Top row: raw data. Bottom row: the same data blurred by the convo- 
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to predict the detailed spectrum of the disk, and the correspondence be- 
tween those two sets of data is by no means clear. However, the models 
can yield data much more closely related to line emissivity than simple 
physical parameters or their combinations. 

To obtain such data, we assume the line emission to originate mainly 
due to the irradiation of the surface layer of the disk by the central white 
dwarf and/or boundary layer (c/. Robinson et al. 1993, Smakl991). Since 
our models are not detailed enough to resolve the surface layer, we assume 
that the lower boundary of the layer coincides with a constant density 
surface Si at which p = pi = 10~®. Typical distance between Si and the 
midplane of the disk, hi was such, that hi/r — 0.08 and hi/r = 0.2 in the 
inner disk and in the outer ring, respectively. Further, we assume that 
the line flux from each element of the layer is proportional to the mass 
contained within that element, pdV , multiplied by the irradiating flux. 
For simplicity, we also assume that all irradiating photons are emitted 
from a point source located at the centre of the white dwarf, so that the 
irradiating flux is given by C/r'^ , where r is the radial coordinate of the 
volume element AV , and C — £„d -I- /Cbi is the combined luminosity of the 
white dwarf and the boundary layer. The distribution of brightness on 
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Figure 4: Tomograms obtained from models Bl and B2 within the irradiation 
approach (blurred in the same way as described in Fig.EJ- 
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Figure 5: Distribution of the compressional power per unit mass averaged over 
the azimuthal angle, — J^^ min ^^^,0^ d^, in models Bl and B2. 



the (vx , Vy ) plane is now given by the integral 

Lirr{Vx,Vy) ^ J p^S{r,9,(j))B{Ux,Vj;,5^^)B{Uy,Vy,5vy)dV. (6) 

The function S{r,0, (f>) describes the shadow cast by Si, and it is given by 

5(r, e,^) = i-n (^j^^ n{p{r\ e, 0) - pi)dr'^ , (7) 

where 

1^ otherwise 

is the Heaviside step function. The third velocity component, Vz, was 
neglected in @ because nearly everywhere in the disk its value was smaller 
than ~ 5% of the local azimuthal velocity, v^. Formally, the integral in 
Eq.Elsubtends the whole space above 5*1. Practically, due to steeply falling 
density, only volume elements closest to 5*1 contribute to it significantly. 
The boundary layer above Si has a mass of ~ 0.15Afdisk, but only about 
25% of its volume is directly illuminated. 

Location and relative intensity of the brightest areas on tomograms 
resulting from the irradiation approach (Fig. |1J agree rather well with 
those observed by Groot (2001) at an advanced outburst phase of U Gem 
(his Fig. 2, Episode 2). The major discrepancy is the bright ring visible in 
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our tomograms at + — 1000 km ^ . In both models the ring is too 
bright compared to the observational data, but relative to the maximum 
intensity obtained in the model it is weaker in B2 where the disk extends 
down to r = 0.045. We conclude that in Bl the ring is enhanced by a 
spurious contribution from the inner edge of the disk at r = 0.1. The ring 
in B2 would be still weaker if absorption of the irradiating flux by the gas 
in the surface layer was taken into account in equation (|HJ . Unfortunately, 
the present models are too crude for such an operation to be reliable. It is 
clear, however, that the effect of absorption should be particularly strong 
for r < 0.1, where h/r is nearly constant (see Fig. |^, and the irradiating 
quanta propagate nearly parallel to Si. Further reduction in ring intensity 
could probably be achieved if the inner boundary of the grid was moved 
even closer to the white dwarf, as the shadow cast by Si at r < 0.045 might 
partly screen the region at r 0.1 where the ring originates. On the other 
hand, in some phases of the activity cycle the intensity of the brightest 
areas in our irradiation tomograms is underestimated relative to the ring. 
This is because substantial pdV work is done by tidal forces directly on the 
gas in the surface layer of the outer disk, where the low-velocity emission 
originates. 

In fact, the heating rate per unit mass, pdV/p, reaches a clear maxi- 
mum just below the surface of the outer disk (see Fig. |3 . For the case of 
U Gem the height of this maximum is ~ 10^^ ergg"^ s~^. With a typical 
outburst accretion rate of M = 3 • 10^* gs~^, and a white dwarf radius 
7?„d = 4-10® cm we get L^i ^ C = \GMiM/R^a ~ 6 • 10^'^ ergs"^ The 
illuminated mass in the region of maximum pAV/ p between r ^ 0.15a 
and r ~ 0.25a approaches ~ 0.025Mdisk, and it is distributed within a 
solid angle of ~ O.Stt. Assuming that the whole incident flux is absorbed 
there we obtain a radiative energy input of ~ 2 • 10^^ ergg^^ s^^, and we 
see that during the outburst the contribution to the line flux from pdV 
heating is rather small. However, pAV dominates just before the outburst, 
when accretion rate is factor of ~ 100 lower. 

The subsurface maximum of pAV/ p in Fig. |3 originates mainly due to 
tidal forcing in the direction perpendicular to the orbit. It also contains 
a contribution from the ambient gas falling onto the disk; however the 
"rainfall" heating is much less efflcient than the tidal one. We checked this 
by re-running simulation Bl with pmin reduced to 3- 10"^*: at f ~ 8.0 Porb 
virtually no changes were seen in the distribution of pAV/ p. 

4 Discussion 

As discussed in Sectional we flnd that spiral waves efficiently propagate 
from excitation regions at the outer edge of the disk toward the white 
dwarf, reaching to at least ~ 0.05a. This conclusion concerns both two- 
dimensional and three-dimensional models; it is however limited to hot 
polytropic disks presented in this paper. The main spiral features seen in 
the density distribution (two-dimensional case) and in the surface density 
distribution (three-dimensional case) are hard to distinguish. On the other 
hand, clear differences are visible in the distributions of the tidal heating 
rate, pAV: in 3D the two main spiral arms are less tightly wound than in 
2D, and a weaker third arm is excited. We suggest that the third arm may 
originate from tidal forcing in the direction perpendicular to the orbital 
plane. The effects of this forcing seem to be responsible for the origin of 
the clear maximum of tidal heating rate per unit mass, pAV/ p, which is 
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located away from the midplane in subsurface layers of the outer disk. 

Doppler tomograms of tidal heating rate derived from both 2D and 3D 
models correlate rather poorly with observed tomograms of U Gem (Groot 
2001). A better agreement (but still not entirely satisfactory) is obtained 
for tomograms of the irradiation flux from the white dwarf through the 
surface layer of the disk. The brightest areas of such tomograms coincide 
with arches observed in U Gem at an advanced stage of the outburst. 
The irradiation tomograms can be derived from 3D models only, which 
indicates that fully three-dimensional modelling is needed for a reliable 
interpretation of the observational data on DN disks. 

According to our results the arches originate in the outer part of the 
disk, fairly high above the midplane {h > O.lr). For this to happen, 
the outer disk would have to be substantially bulged. The bulging phe- 
nomenon may be explained within the following (not entirely new) sce- 
nario: prior to the outburst the gas transferred from the secondary mainly 
collects in a ring at the circularization radius, and only partly accretes 
through the disk onto the white dwarf. The ring expands as the gas flows 
in, but it remains cool until heating from tidal forcing in the orbital plane 
grows so strong that it cannot be balanced by radiative cooling. The 
ring begins to expand even more rapidly, and within it the gas located 
away from the midplane begins to receive additional internal energy from 
tidal forcing in the direction perpendicular to the orbital plane. Eventu- 
ally, strong spiral shocks develop, and a dynamical instability of the kind 
described by Rozyczka & Spruit (1993) sets in. 

Obviously, such scenario cannot be the whole story, as it is not linked 
to the thermal instability believed to be at least partly responsible for 
eruptive phenomena in DN and other classes of cataclysmic variables. 
Nevertheless, it seems to indicate a promising direction of further research. 
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